function [R,E]=mcf(A,epsilon)
% Modified Chelosky factorization of indefinite matrix
% Based on Algorithm MC in Practical optimization (GMW 81) p111
% INPUT:
%       A - the matrix to factorize;
%       epsilon - a small positive numeber (machine precision);
% OUPUT:
%       R - the triangular factor;
%       E - the error matrix R'R=A+E;
% REVISIOn:
%       jqin, 01/feb/2011, created;
gamma=max(diag(A));
xi=max(max(A-diag(diag(A))));
n=numel(A(:,1));
nu=max([1,sqrt(n^2-1)]);
beta=sqrt(max([gamma,xi/nu,epsilon]));
R=spalloc(n,n,(1+n)*n/2);
E=spalloc(n,n,n);
mu=zeros(n-1,1);
for k=1:n
    if k~=n
        mu(k)=max(abs(A(k,k+1:n)));
        R(k,k)=max([epsilon,sqrt(abs(A(k,k))),mu(k)/beta]);
        E(k,k)=R(k,k)^2-A(k,k);
        for j=k+1:n
            R(k,j)=A(k,j)/R(k,k);
            for i=k+1:j
                A(i,j)=A(i,j)-R(k,j)*R(k,i);
            end
        end
    else
        R(k,k)=max([epsilon,sqrt(abs(A(k,k)))]);
        E(k,k)=R(k,k)^2-A(k,k);
    end
end